
program define rdbinscatter
	version 14
	syntax varlist(min=2 max=2) [if] [, CUTpoint(real 0) binsize(integer 10) ///
		xname(string) yname(string) func(string) raw_data(integer 1) ///
		start(real 0) end(real 0) tick_dist(real 0)]
	local y = word("`varlist'", 1)
	local x = word("`varlist'", 2)
	if "`func'"=="" local func = "lfit"
	preserve
	cap keep `if'
	cap keep if `y'!=. & `x'!=.
	gen __xx__rv = `x' - `cutpoint'
	gen __xx__abs_rv = abs(__xx__rv)
	gen __xx__treat = __xx__rv>=0
	gsort __xx__treat __xx__abs_rv
	by __xx__treat: gen __xx__bin = ceil(_n/`binsize')
	bys __xx__treat __xx__bin: gen __xx__obs = _N
	egen __xx__mean_x = mean(`x'), by(__xx__treat __xx__bin)
	egen __xx__mean_y = mean(`y'), by(__xx__treat __xx__bin)
	if `start'==0 & `end'==0 & `tick_dist'==0 {
		if `raw_data'==0 twoway (scatter __xx__mean_y __xx__mean_x if __xx__obs==`binsize', mc(gs2)) ///
			(`func' `y' `x' if __xx__treat==0, lc(gs2)) ///
			(`func' `y' `x' if __xx__treat==1, lc(gs2)), ///
			xti("`xname'") yti("`yname'") ///
			legend(off) xli(`cutpoint', lc(gs10) lp(dash)) ///
			graphregion(color(white))
		if `raw_data' twoway (scatter `y' `x', mc(gs10) msi(vsmall)) ///
			(scatter __xx__mean_y __xx__mean_x if __xx__obs==`binsize', mc(gs2)) ///
			(`func' `y' `x' if __xx__treat==0, lc(gs2) alw(none) fc(gs10%25)) ///
			(`func' `y' `x' if __xx__treat==1, lc(gs2) alw(none) fc(gs10%25)), ///
			xti("`xname'") yti("`yname'") ///
			legend(off) xli(`cutpoint', lc(gs10) lp(dash)) ///
			graphregion(color(white))
	}
	else {
		if `raw_data'==0 twoway (scatter __xx__mean_y __xx__mean_x if __xx__obs==`binsize', mc(gs2)) ///
			(`func' `y' `x' if __xx__treat==0, lc(gs2)) ///
			(`func' `y' `x' if __xx__treat==1, lc(gs2)), ///
			xti("`xname'") yti("`yname'") ///
			legend(off) xli(`cutpoint', lc(gs10) lp(dash)) ///
			xlab(`start'(`tick_dist')`end') xsc(r(`start' `end')) ///
			graphregion(color(white))
		if `raw_data' twoway (scatter `y' `x', mc(gs10) msi(vsmall)) ///
			(scatter __xx__mean_y __xx__mean_x if __xx__obs==`binsize', mc(gs2)) ///
			(`func' `y' `x' if __xx__treat==0, lc(gs2)) ///
			(`func' `y' `x' if __xx__treat==1, lc(gs2)), ///
			xti("`xname'") yti("`yname'") ///
			legend(off) xli(`cutpoint', lc(gs10) lp(dash)) ///
			xlab(`start'(`tick_dist')`end') xsc(r(`start' `end')) ///
			graphregion(color(white))
	}
	restore
end
